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Standard methods of structural dynamic analysis assume that the structural characteristics are deterministic. 
Recognizing that these characteristics are actually statistical in nature, researchers have recently developed a 
variety of methods that use this information to determine probabilities of a desired response characteristic, such 
as natural frequency, without using expensive Monte Carlo simulations. One of the problems in these methods 
is correctly identifying the statistical properties of primitive variables such as geometry, stiffness, and mass. 
We present a method where the measured dynamic properties of substructures are used instead as the random 
variables. The residual flexibility method of component mode synthesis is combined with the probabilistic methods 
to determine the cumulative distribution function of the system eigenvalues. A simple cantilever beam test problem 
is presented that illustrates the theory. 


Nomenclature 

b = boundary degrees of freedom 

[C] a = correlation matrix of random variables in 

substructure a 

g = limit state function 

i = internal degrees of freedom 

k = number of kept modes per substructure 

[L]“ = Cholesky decomposition of [C] a 

N = total number of degrees of freedom 

per substructure 

p = total number of substructures 

Pf = probability of failure 

p(x < X) = probability that x is less than X 

{«} = vector of uncorrelated standard normal random 

variables 

X - vector of random variables 

X * = design or most probable point 

(x) = vector of correlated normal random variables 

{x}' = vector of correlated standard normal random 

variables 

a - substructure number 

Introduction 

S TRUCTURAL analysts have always known that the parameters 
of the system being modeled are not deterministic, because of 
manufacturing tolerances, material deviation, and other factors. Un- 
til recently, the primary way to deal with this knowledge was to use 
safety factors, which are qualitative and based primarily on experi- 
ence. In an effort to take account of these variations in the structural 
parameters in a more quantitative fashion, significant research has 
been performed to develop methods to actually use the statistical 
characteristics of the input quantities in the analysis to generate an 
output value that is also described statistically. Monte Carlo (MC) 
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simulations can be performed to calculate these probability distri- 
butions, but up to a million runs are required for accurate results. 
Approximate techniques have therefore been developed that require 
several orders of magnitude fewer calculations than MC techniques. 
One such method, the fast probability integration (FPI) method, 1 
has recently been implemented in a new probabilistic finite element 
code, NESSUS. 2 

Numerical analyses of structural vibration generally use the finite 
element method (FEM) as the basis for obtaining free and forced 
response characteristics. A frequent problem when using the FEM 
is that, for large models composed of many substructures, the num- 
ber of degrees of freedom (DOFs) is so large that the computational 
costs are prohibitively expensive for eigenanalysis. This problem is 
particularly relevent when using probabilistic techniques, because 
the eigenanalysis has to be repeated many times. The favored so- 
lution method in industry for the deterministic situation is to apply 
dynamic component mode synthesis (CMS) methods, which sub- 
stantially reduce the DOFs for the system model. 

This paper defines a procedure for combining CMS with proba- 
bilistic methods to obtain the statistical characteristics in an efficient 
manner. These characteristic are summarized in the form of the cu- 
mulative distribution function (CDF). The procedure makes use of 
statistical distribution information of each substructure’s dynamical 
modes and residual flexibility, which are available from experimen- 
tal data. This information is synthesized into a system model using 
the residual flexibility method of CMS, and the statistics of the sys- 
tem dynamic characteristics are obtained using FPI. An advantage 
of this method over existing probabilistic structural analysis meth- 
ods is that, in many cases, the statistics of the substructure dynamic 
characteristics may be easier to determine than those of primitive 
random variables (RVs) such as geometry, material stiffness, or den- 
sity. In addition, these statistical dynamic characteristics completely 
incorporate all random factors in the structure. This is virtually im- 
possible to do with primitive RVs such as thickness or density, since 
they are in reality random fields, varying over the entire structure. 
Final development of the method should allow probabilistic meth- 
ods to be applied to much larger models than previously possible, 
such as turbomachinery bladed disks, which are composed of many 
almost identical substructures whose structural characteristics can 
be described statistically. 

Probabilistic Theoretical Background 

Research in the field of probabilistic structural mechanics has 
concentrated in two areas. The first can be described as perturbation 


830 



BROWN AND FERRI 


831 


methods, and the second as reliability methods. The perturba- 
tion method, as developed by Collins and Thomson, Kiefling, and 
Collins, Kennedy, and Hart, is used to derive an analytical expres- 
sion for the mean and standard deviation of structural eigenvalues 
and eigenvectors as a function of the derivatives of the mass and 
stiffness matrix for each input RV. 3-5 Hasselman and Hart 6 used 
CMS to derive analytical expressions for the system, or global, 
eigenvalue matrix as a function of the modally reduced substruc- 
ture stiffness and mass matrices. This approach is presently being 
studied by Mahadevan and Mehta. 7 

The research described in this paper employs the reliability 
method approach to determining the statistical structural response 
characteristics. To review this technique briefly, consider the scalar 
limit state function g(X), defined as some value of interest that is 
a function of a vector of RVs. Cornell, 8 using the reliability ap- 
proach, developed what would later be called the first-order relia- 
bility method (FORM) by truncating the higher-order terms (HOTs) 
from the series 



of g, Hasover and Lind introduced an initial reduction of each of 
the primitive normal RVs Xj to standard normal RVs u f using 


= g(v-x) + — l*i) + HOTs (1) 

7Z 1 !lx ‘ 

where is a vector of the mean values of the RVs in X. This 
resulted in the following simple approximations for the mean value 
and standard deviation (cr) of g: 


M* ^ g(t*x) 



( 2 ) 

( 3 ) 


Hasofer and Lind 9 refined and expanded this method further. They 
redefined g as 


g(X) = Y(X)-y 


(4) 


where y is a specific value and the performance function Y(X) is 
a function of the RVs. This formulation divides the space into two 
parts, g < 0(Y < y) and g > 0 (Y > y). The probability that the 
function Y does not exceed the value y is the probability that g < 0. 
For example, if y equals an eigenvalue of interest, p(g < 0) will 
be the probability that the actual eigenvalue obtained is less than 
the one of interest. There will therefore be a limit state function 
g for every y in the range of possibility. If the distributions of the 
RVs are normal and the resulting distribution of g is normal, then 
the following steps can be performed to obtain the probability that 
g < 0 (Y < y). First, a transformation of g to standard normal 
coordinates Z is made as follows: 


g(X) - 


(5) 


and if p, which Cornell termed the safety index, is defined as 


This transformation allows the new limit states g(V) = 0 to 
be plotted in standard normal joint probability space for every 
possible y. The joint probability of g < 0 will be the volume un- 
der the multidimensional bell-shaped surface over the area where 
g < 0. Furthermore, p can be shown to be the shortest distance 
from the origin to the line g(V) = 0, and the point on the line at p 
can be called the most probable point (MPP), or design point X\ 
because that point will have the highest probability of occurrence 
of any point along the line g(V) = 0 (see Fig. 1). 

The reliability method was expanded by Rackwitz 10 to multi- 
dimensional, nonlinear limit states in which g is an explicit func- 
don of the RVs. Wu and Wirshing 1 1 developed the advanced mean 
value (AMV) method, a procedure for minimizing the number of 
iterations needed to obtain the MPP for nonlinear limit states; this 
is vital for nonexplicit limit states, such as finite element solutions. 
The limit state is approximated as a linear function about the means 
of the RVs [Eq. (1)], and the partial derivatives are approximately 
obtained by numerically differentiating the limit state with respect 
to each RV. Values of p and X* are obtained for each desired limit 
state by using the FPI method, which is a compilation of the im- 
provements to the FORM made by Wu, 11 Rackwitz, 10 and others. 
At this point in the procedure, an exact solution for each of the 
limit states (usually a finite element solution) is found by plugging 
in these most probable points. These results and their associated p 
values can then be used to create an entire CDF, which is shown in 
Wu’s paper 11 to be in very good agreement with MC simulations 
for several examples. Further iterations can be performed by ex- 
panding the limit state about the new design points, instead of about 
the means as in Eq. (1). The FPI and AMV methods have been 
incorporated in NESSUS, 2 a probabilistic finite element program 
under development by NASA Lewis and the Southwest Research 
Institute. 


P = V-s/Vg ( 6 ) 

then 

p[g(X) < 0] = p(Z < -P) = <t>(-p) (7) 

where <J> is the CDF of the standard normal distribution function 
found in handbooks. Since negative values for p are not tabulated, 
the relationship 

<t> ( _0) = !_<*>(£) (8) 

is used instead to calculate this probability. The complete CDF is 
formed by finding the p values for all the limit states in the range of 
possibility. These values are only accurate, though, for linear limit 
states where the RVs have normal distributions. In addition, g could 
be formulated differently for some cases, thereby yielding different 
probabilities. To provide invariance with respect to the formulation 


Probabilistic Dynamic Synthesis 

The proposed methodology makes use of the residual flexibility 
method of CMS. This method has been developed by MacNeil, 12 
Craig and Chang, 13 and Martinez et al. 14 The essential idea in CMS 
is that substructure modes are truncated because their higher modes 
will not have a major effect on the system modes. The residual 
flexibility method incorporates the effects of the higher modes by 
determining their flexibility. A side benefit is that all the elements 
of the system stiffness matrix can be obtained from test and that the 
mass matrix can be closely approximated by a unity matrix in the 
nonboundary partition (equal in size to the number of kept modes, 
k ). Since all the substructure information can be obtained from test, 
probabilistic data can be completely incorporated into the system 
matrices to obtain the system modes. 

The first step of the probabilistic dynamic synthesis (PDS) method 
developed in this paper is to divide the model of a structure into sub- 
structures or = a, b , . . . , p. The physical displacement vectors of 
each substructure, which have either a subscript i denoting internal 
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DOFs or a subscript b denoting boundary DOFs, can be written 
as 



where there are altogether dim*,- + dim** = N DOFs for that 
substructure. 

Each substructure is represented by n samples. Each sample 
is modally tested individually in a configuration such that the 
interface locations with other substructures are in a free condi- 
tion. For substructure a, sample i, the test will yield eigenvalues 
A® and eigenvectors {^?“}. In addition, the boundary partition of 
the residual flexibility matrix [ G is obtained from the measured 
boundary drive- point frequency response functions of the boundary 
coordinates. 15 For use in the FDS method, only the kept (nontrun- 
cated) eigenvalues, the boundary coordinates of the kept eigenvec- 
tors, and the boundary partition of the residual flexibility matrix are 
needed. These values can be combined into a single vector {x}?, 
defined as 

tM'V 


{x}? = 


x' 


(ii) 


of [L]“ and then into the original RVs {x} using Eq. (9). The new 
matrices {A}, [<J>], and [G bb \ are pulled out from (x) and placed in 
substructure mass and stiffness matrices according to the residual 
flexibility formulation: 


K a = 


A + <t-LGj<t> bk -<t>I k c;,y 


^bk^bb ' 

sym 


J bb 



where A is a diagonal matrix of the eigenvalues {A}. 

The system mass and stiffness matrices are now generated by 
directly coupling the substructure mass and stiffness matrices. This 
is accomplished by ordering the kept DOFs of each substructure 
sequentially in the system matrices and adding the boundary par- 
titions together. The system eigenvalues are then obtained, and a 
single eigenvalue of interest is chosen. As each independent RV in 
the set of p vectors {«} is varied, the numerical partial differentiation 
of the eigenvalue with respect to the RV can be calculated. These 
are used to generate the linear approximation of the limit state at 
the mean value [Eq. (1)], which is used by the FPI code to obtain 
the MPPs and a first mean- value solution for the CDF of the chosen 
system eigenvalue. The MPPs are plugged back into the substruc- 
tures’ mass and stiffness matrices, the system is resynthesized, and 
new, updated eigenvalue levels are obtained for each probability 
level, following the AMV method described by Wu and Wirshing. 11 
These levels are then plotted to show the entire CDF. 


*• 

l J , 

where {ipi,} j is a vector of the boundary-node modal displacements 
for the yth mode, and X J is the ;th eigenvalue of substructure or. 

If the entire sample of substructure or is tested, {x} a can there- 
fore be defined as a vector composed of elements that are each a 
normally distributed RV with measured mean and standard devi- 
ation. Using Eq. (9), this vector is now converted to {*)'“, a vec- 
tor of standard normally distributed RVs. In addition, there will 
be some degree of correlation between the RVs. These correlation 
values range from 0, or no correlation, to ±1, or fully correlated, 
and can be easily calculated from the measured data. The values 
are placed in a correlation matrix [C] a relating each element with 
every other element. For the probabilistic analysis, a set of inde- 
pendent standard normal RVs {w}° will have to be obtained. This 
can be accomplished 16 by making an orthogonal transformation of 
[x} ,a with its Cholesky-decomposition lower triangular matrix [L]“ 
to uncouple the {x}' coordinates, thereby creating {w}®. This can be 
expressed for substructures a = a, b, . . . , p as 

{*r = uw a 02) 


Test Case 

Analysis of a spring-mass system (Fig. 2) using the PDS method 
has been completed. The test system consists of two substructures, 
a and b , each having four DOFs. Five thousand samples of each 
substructure were created initially using standard MC techniques. 
To achieve complete probabilistic generality, each spring in the sys- 
tem was assigned a normal distribution with a mean of 200 and 
standard deviation of 10, and each mass was assigned a normal dis- 
tribution with a mean of 1.0 and standard deviation of 0.5. Since 
this initial numerical simulation is performed to represent measure- 
ments of the dynamic characteristics of a physical population, any 
distribution can be chosen for the masses and springs as long as it is 
consistent for the full-up model MC case and the PDS case. Normal 
distributions were therefore chosen for convenience. The effect of 
nonnormal distributions will be discussed later. The MC random 
vectors were then used to create the mass and stiffness matrices for 
the substructures (5000 for each), and a modal analysis run on the 
substructure samples to obtain their eigenvalues {A}“ and eigen- 
vectors [<t>] a . Three of the four modes for each substructure were 
kept for the analysis. The boundary partition of the N x N residual 
flexibility matrix [G h b] a was analytically calculated directly from 
the modes that had been chosen to be truncated, in this case just the 
highest one, using 


The FPI algorithm requires that each independent RV be varied 
individually by some percentage of its standard deviation a , which 
was chosen to be 50% for this development, while the other RVs 
are kept constant at their mean values. The choice of the size of the 
variation is somewhat arbitrary. Each of these cases is then back- 
transformed to form a corresponding case of the original correlated 
RVs. These are then plugged into the model to generate the limit 
state approximation [Eq. (1)] of the response value, which is used 
to obtain the $ values and the design points X * . 

Since the distributions of the RVs are standard normal, 0.5a will 
simply equal 0.5 for the RV to be varied. The first case is therefore 


[C]= ^ MM (is) 

i — k + I * 

where, in this case, k = 3 and N = 4. The statistics on these dy- 
namic characteristics and the correlation between them were then 
calculated. These statistics are listed in Table 1 . The listed quantities 
comprise the vectors {x}“ and {x} /; as described in Eq. (1 1). 

A distribution characterization routine 17 was also performed on 
the distributions to see if they could be characterized as normal, 
which is an assumption of the methodology outlined. Partial results 
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The next case will consist of the second element in {«}“ equaling 
0.5 and all the other elements of [u}° as well as all the elements of 
the other (u) m equaling zero, and so on. 

For each case, the {«} for each substructure is then transformed 
to the set of correlated standard normal RVs (x)' using the transpose 
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Fig. 2 Test-case system. 
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Table 1 Statistics of dynamic characteristics 


Eigenvalues Mean 

Standard deviation 


Substructure A 


1 

30.431 

1.2702 

2 

246.74 

10.336 

3 

552.89 

23.014 


Eigenvectors, boundary location only 

1 

0.7073 

0.01155 

2 

-0.70725 

0.0285446 

3 

0.70688 

0.058827 


Residual flexibility ( one boundary point only) 

l 

6.46E-04 

1.6159E— 04 


Substructure B 


1 

0.0 

0.0 

2 

150.51 

6.706 

3 

489.03 

21.916 


Eigenvectors, boundary location only 

1 

0.53466 

6.9105E— 03 

2 

0.75592 

2.1945E— 02 

3 

-0.55443 

7.4042E— 02 


Residual flexibility (one boundary point only ) 

1 

6.46E-04 

1.5347E— 04 


Table 2 Partial results of distribution-types 
routine: substructure B y eigenvalue 2 


W statistics (all types included ) 

Normal 

0.00948 

Exponential 

0.30339 

Weibull 

0.04448 

EVD 

0.04323 

Lognormal 

0.00934 

Normal distribution parameters 

Sample mean = 

150.61 

Sample std. dev. = 

6.502 

Normal-distribution CDF fit to data 

Response value 

CDF value 

143.0 

0.1210 

147.0 

0.2895 

150.0 

0.4628 

155.0 

0.7503 

Lognormal parameters, base e 

\x = 5.014 


a = 0.43201 

Lognormal-distribution CDF fit to data 

Response value 

CDF value 

143.0 

0.1193 

147.0 

0.2947 

150.0 

0.4713 

155.0 

0.7540 


of the distribution characterization routine for one of the RVs are 
shown in Table 2. The W statistic is a goodness-of-fit test developed 
by Wirshing and Carlson, where a smaller number indicates a closer 
fit. lx The results show that the data are well represented by both 
normal and lognormal distributions, with the lognormal being a 
little better. The CDF values for each distribution, however, indicate 
that the curves for the two distributions predict almost exactly the 
same value, which can be the case for a particular set of lognormal 
parameters. The assumption of normality was therefore deemed to 
be accurate. 

At this point, the procedure follows the outline discussed previ- 
ously for an actual case, which would use modal testing of physi- 
cal samples to generate the dynamic characteristics and correlation 
matrix. For each substructure, a matrix composed of cases of {«} 
vectors [see Eq. (14)] was generated and multiplied by the correla- 
tion Cholesky decomposition matrix [L] c to obtain {*}', the set of 
correlated standard normal RVs. These were then converted to their 
nonstandard normal distributions and used in the residual flexibility 
substructure stiffness matrix. The substructures were then coupled 


Table 3 Sample MPP output 

Fund, eigenvalue response value = 8.101 
Probability = 0.01 
Most probable point 
Substructure-a RV Value 



-0.4863 

u 2 

-0.1692 

«3 

-0.05685 

«4 

-1.5807 

"5 

-0.1221 

U6 

-0.0269 

«7 

0.211 

Substructure-^ RV 

Value 

Ml 

-1.5205 

U 2 

0.2641 

M3 

-0.2361 

M4 

-0.0085 

"5 

-0.3057 


-0.0473 

«7 

0.2267 


Table 4 Dispersion error of PDS vs full MC model 


CDF 

value 

AMV 

A from 
AMV mean 

MC 

A from 
MC mean 

Error 

of A, % 

0.010 

8.116 

-0.635 

8.094 

-0.613 

-3.46 

0.050 

8.300 

-0.451 

8.266 

-0.441 

-2.23 

0.100 

8.399 

-0.352 

8.352 

-0.355 

0.87 

0.159 

8.476 

-0.275 

8.426 

-0.281 

2.22 

0.200 

8.520 

-0.232 

8.467 

-0.240 

3.62 

0.500 

8.751 

0.000 

8.707 

0.000 

Not applicable 

0.800 

8.985 

0.234 

8.944 

0.237 

1.44 

0.841 

9.029 

0.278 

8.987 

0.280 

0.74 

0.900 

9.108 

0.357 

9.072 

0.365 

2.25 

0.950 

9.210 

0.459 

9.169 

0.462 

0.65 

0.990 

9.403 

0.652 

9.344 

0.637 

-2.28 



Fig. 3 CDFs for fundamental eigenvalue using PDS method and MC 
of full, unsubstructured model: — , pds amv and , full mod. me. 


together, and a modal analysis was performed on the system ma- 
trices. The first system eigenvalue for each case, which was the 
response value chosen, was then input along with its {«} case into 
the FPI algorithm routine. 

The output of the FPI routine is the mean- value solution, an initial 
estimate of the CDF of the response variable, and the MPPs for 
the specified CDF probability levels. One MPP from the output 
is shown in Table 3. Following the AMV procedure, these MPPs 
were recorrelated and converted to the original dynamic RVs as 
before, coupled, and a solution obtained for the updated fundamental 
eigenvalue. This value and its associated probability level were then 
used to create a new CDF. 
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For verification of the PDS method, an MC analysis was per- 
formed on the same nondeterministic spring-mass system, with the 
system eigenvalue directly obtained from each unsubstructured sam- 
ple. The CDF for this full model is superimposed on the AMV CDF 
from the PDS method in Fig. 3. A very small amount of error is 
indicated graphically. To identify the error quantitatively, the amount 
of variation of the fundamental eigenvalue from its mean value at 
selected probability values for the PDS method was compared with 
the spread for the full model. The result, shown in Table 4, indi- 
cates that the deviations from the mean as computed by the AMV 
and MC methods agree to within 5%. In addition, the mean value 
of the fundamental eigenvalue computed by the AMV method is 
8.75 1 , which is only 0.5% higher than the one computed using MC 
(8.707), and the AMV standard deviation is 0.276, which is only 
1 .8% less than the MC standard deviation (0.281). 

Concluding Remarks 

A new methodology has been presented for performing analysis 
of structures composed of substructures whose dynamic characteris- 
tics can be identified statistically. This method uses the substructure 
eigenvalues, eigenvectors, and residual flexibility as random vec- 
tors for determining the desired response value by combining new 
probabilistic analysis techniques with the residual flexibility method 
of CMS. By this combination, improved computational efficiency 
is obtained as well as more confidence in the input statistics. Re- 
sults for a test case show the method predicts close to the same 
CDF for the system fundamental eigenvalue as a non substructured 
probabilistic MC analysis. 

Future work on this method will include examining some ba- 
sic conceptual questions on the limitations and applicability of the 
method. The fact that the dynamic RVs exhibited normal distribu- 
tions may be a result of the assumption of normal distributions for 
the original primitive RVs. One extension of this research will be to 
examine the effect of nonnormal primitive RVs, which could result in 
nonnormal dynamic RVs. Equivalent normal distributions for these 
RVs would have to be obtained using methods developed by Chen 
and Lind and other researchers. 19 The methodology for handling the 
correlation of these nonnormal RVs will also have to be examined, 
as the correlation coefficient may not adequately describe the rela- 
tionship between the variables. Since the research described in this 
paper is the first step in developing an entirely new methodology, the 
assumption of normal primitive RVs was necessary to avoid these 
complexities, and since normal distributions are commonly found, 
the results are applicable and worthwhile. 

Another topic of interest is the formulation of a hybrid method 
combining analysis and test, since some of the dynamic charac- 
teristics necessary for the synthesis may be difficult to measure in 
some circumstances. Other questions include finding cases where 
the number of RV s can be reduced by perhaps only allowing stiffness 
or mass to vary, and examining the effect of boundary variability for 
situations like the fir-tree interface between blades and disks, which 
are neither fixed nor free. In addition, the method will be compared 


with perturbation methods to determine the areas of most efficient 
applicability for each one. 
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